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Abstract 

We present an efficient numerical technique to evaluate the matrix of the (quasiparticle)-random- 
phase approximation, using the finite amplitude method (FAM). The method is tested in calculation 
of monopole excitations in ^^'^Sn, compared with result obtained with the former iterative FAM. 
The neutron-pair-transfer modes are calculated with the present method and their character change 
in neutron-rich Pb isotopes is discussed. Computational aspects of different FAM approaches are 
also discussed for future applications to a large-scale computation. 

PACS numbers: 21.60.Jz ; 21.10.Re ; 24.30.Cz 
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I. INTRODUCTION 



The random-phase approximation (RPA) is a leading theory in studies of elementary 
excitations in nuclei and other quantum many-body systems. The RPA, which is equivalent 
to the small-amplitude limit of the time-dependent mean-field theory, is commonly known 
in the matrix form jl, 2| as 

-B* -A* I \ / " \ 



(1) 



Here, X" and are respectively called forward and backward amplitudes of the n-th RPA 
normal modes. We refer to the matrix in the left-hand side of Eq. ([1]) as the RPA matrix, or 
as the quasiparticle RPA (QRPA) matrix in the case that the mean fields contain the pair 
potential. In nuclear mean-field models, the evaluation of the (Q)RPA matrix is a tedious 
task and requires significant efforts for programing the computer code. The main purpose 
of the present paper is to present a feasible and efficient numerical method to evaluate the 
(Q)RPA matrix. 

The finite amplitude method (FAM) was proposed as a feasible numerical approach to the 
calculation of the strength functions [3|] . Recently, it has been extended to the quasiparticle- 
basis representation to include the pairing correlations jj]. In this method, instead of cal- 
culating the RPA normal modes, one solves the linear response equation with an external 
field at a given frequency. The FAM has been successfully applied to the Skyrme energy 
functionals in different representations; the three-dimensional (3D) coordinate-mesh repre- 
sentation with no symmetry restriction but no pairing {5,6], the quasiparticle basis with the 
radial-mesh representation for spherical nuclei and that with the harmonic-oscillator- 



basis representation for axially deformed nuclei |7||. The FAM is a 



so shown to be superior 
7|. 



to the conventional approach, with respect to numerical costs 

So far, the FAM has been utilized for the calculation of the strength functions, using 
iterative algorithms. Hereafter, this approach is referred to as "iterative FAM" (i-FAM). 
However, when we are interested in low-lying discrete modes of excitation, it is desirable to 
obtain the RPA normal modes in Eq. ([T]). In this paper, we show another usage of the FAM, 
for explicit construction of the RPA matrix. We call this approach "matrix FAM" (m-FAM) 
in this paper. The method only requires a straightforward extension of the former FAM for 
the strength function. Since the complicated programing is required for calculation of the 



RPA matrix for realistic energy functionals, this new method is useful for the verification of 
the existing/developing computer codes of the (Q)RPA, as well. 

The paper is organized as follows: In Sec. HIl after a brief introduction to the finite 
amplitude method, we propose a new method, m-FAM, of calculating the RPA matrix 
for systems without pairing correlations. In Sec. IIIH we present a method applicable to 
superfiuid systems, namely the m-FAM for the QRPA. Numerical results are shown in Sec. IIVI 
for monopole strength function in ^^''Sn and the neutron-pair transfer to excited 0^ states 
in neutron-rich Pb isotopes. Comparisons between i-FAM and m-FAM are shown in Sec. |Vl 
in terms of the computational point of view. Finally, the summary is given in Sec. IVI[ 

II. FAM CALCULATION OF THE RPA MATRIX 

First, let us discuss the case without the pairing correlation. It is customary to use 
the canonical single-particle representation, which we adopt in this paper too. Thus, the 
forward and backward amphtudes in Eq. ([T]) have particle-hole (ph) indices, {Xph,Yph), 
in which the particle and hole states are assumed to be eigenstates of the single-particle 
Hamiltonian at the ground states: h[pQ]\(f)k) = €k\4>k)- We use the notation (f)h for hole 
orbitals {h = 1, ■ ■ ■ , A) and 0p for particle orbitals {p = A + 1, - ■ ■ , oo). The matrices A and 
B in Eq. ([T]), which have two pairs of ph indices, are given by 

^ph,p'h' — [^p — ^hjOpp'Ohh' + , ^ph,p'h' — ^ • \^) 

Opp'h' OPh'p' 

The residual interactions, dh/dp, are the derivatives of the single-particle Hamiltonian h[p] 
with respect to one-body density, evaluated at the ground-state density (p = po). The 
explicit evaluation of these residual interactions is the most demanding part in the RPA 
calculations, with respect both to the computational cost and to the programming task. 
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A. The Finite Amplitude Method (FAM) 

In this subsection, we recapitulate the FAM js^. The upper part in the left-hand side of 
Eq. leads to 

{Aph,p'h'^p'h' + Bphpih'Ypfhi) = (Cp — th)Xph + (3) 

p'h' 

«'.'-E(l!^Av. + |^y„..). (4) 



p'h' 

f. Isi 



. dpp'h' dph'p' , 

In Ref. (3|, we have proposed the FAM which provides an easy way to evaluate bhph for a 
given vector (X, y). All we need to do is to calculate the single-particle Hamiltonian h\p\ 
at the density slightly different from the ground state p = po as, follows. Provided that a 
real parameter rj is small enough to allow us to neglect the rf and higher-order terms, 

5/1 = i {h\p^] - h[po]) , (5) 

Pri = Po + r]6p = \A){'^h\ + 0{rf), (6) 

h 

where the hole orbitals are slightly modified from the ground-state canonical states in 
different manners between the ket and bra states, as follows: 

\i)h) = \<Ph) +r]YXph\(pp), {'4}h\ = {<ph\+r]YYph{(pp\. (7) 

n 

Note that the use of different bra's and ket's in Eq. (I6l) leads to non-hermitian 5p [3[. From 
Eq. ([6]), one can easily see that this is necessary to obtain the well-known relation in the 
RPA; 5pph = TjXph and 5php = V^ph- 

Now, using the finite-difference calculation of Eq. ([5]), we can evaluate Shph for a given 



3,3, 



the linear response 



(X, Y). In earlier works on calculations of the strength functions 
equation was solved with an iterative algorithm starting from an arbitrary initial vector 
{X, Y). The convergence of the iteration provides the self-consistent RPA amplitudes {X, Y). 
In this iterative process, the RPA matrix elements themselves, A and B, are never calculated, 
instead, only the product of the RPA matrix and the vector, such as AX+BY, are calculated. 

B. Calculation of the RPA matrix 

Now, we present the essential idea of the present paper, that is a method of calculating 
the RPA matrix without explicit evaluation of the derivative, dh/dp. The idea is very simple 
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and immediately understood from Eq. Namely, if we choose the vector as 



-^mi ^mp' ^ih' 1 ^^mi 0, (8) 

where m > A and i < A, then, the FAM calculation of 6hph in Eq. ([5]) leads to dhph/dpp'h'- 
If we choose 

then, it produces dhph/dph'p'- Since the FAM can provide the vector {AX +BY)ph for a given 
{X,Y), we can obtain the matrix elements, Ap^^p'h' and Bph^p'h', by choosing the vectors, 
Eqs. ([8]) and iQ, respectively. Therefore, the RPA matrix can be explicitly constructed by 
the calculation of the single-particle Hamiltonian h[pri] only, using the FAM. 

Sometimes, we resort to a different from of the RPA equation. In such difference 
choice of the vector may be convenient. For instance, the RPA equation ([1]) can be recast 
into the canonical form in terms of the normal- mode coordinate Q*-"^ and momentum P^"^ 

(A + S)(A-5)g(") = w2g(n)^ (A-5)(A + 5)P(") =a;^p("\ (10) 

Mapping these normal coordinates on given collective variables, we may obtain the RPA 
(Thouless-Valatin) collective mass, which has been recently used to study the collective 
quadrupole dynamics 8|-[1Q|]. The numerical solution of the RPA eigenvalue problem can 
be simplified by further transforming Eq. flTUl) into its hermitian form An advantage of 
;his approach in the large-scale parallel computing has been recently demonstrated as well 



These calculations do not require the matrix A and B separately, but need the sum of 



them, {A ± B)ph^p'h'- This is directly accessible with the FAM choosing (X, Y) as follows: 

To construct the full RPA matrix, the residual fields should be calculated, according to 
Eq. (|5]), for all the independent unit vectors, Eqs. ([8]) and Q. However, the present m-FAM 
does not resort to iterative solver. Thus, if the dimension of the RPA matrix is small, the 
m-FAM is computationally more efficient than the i-FAM. The detailed discussion about 
the computational aspects is shown in Sec. El 



5 



III. FAM CALCULATION OF THE QRPA MATRIX 



In this section, the result in Sec. |TT]is generahzed for the QRPA with the pairing corre- 
lations. The mean-field Hamiltonian in the ground state is diagonalized in the quasiparticle 
states: 

H = J2 ( (^'^^ - ^^kiHci + ^Akicl4 + ^AliQCjA =Eo + J2 E^.a\a^. (12) 

kl ^ ^ ii 

Here, (cfc, c|^) are the annihilation and creation operators of a particle at the basis state 
fc, and (a^, aj^) are those of the quasiparticle states. The single-particle (ph) Hamiltonian 
/ifci[p, K, K*\ and the pair (pp,hh) potential Afczfp, k, k*\ are now functionals of the one-body 
density p^/ and the pair tensors {nkhi^li)- "^^^ forward and backward amplitudes have 
two-quasiparticle (2qp) indices, X^,^ and Y^u, and the QRPA matrix has a pair of 2qp 
indices, A^y^^iyi and B^^^^iyi. Here, the quasiparticle states {Uk^, V^^) are chosen to be states 
corresponding to the positive quasiparticle energies E^. 

A detailed formalism of the FAM for the QRPA is found in Ref. j4|. Here, we briefiy 
summarize the main result. Again, the upper part in the left-hand side of Eq. ([T]) is written 



as 



{A^u^fjt'y'X^iyi + B^y^^i^iY^iyi) = {E^ + Ey)Xfj_y + SHj^^\ (13) 
SHll = {U^hV* - V^6A^-'>*V* + f/^(5AW[/* - V^h^U*)^^ , (14) 
where the induced fields, 6hki and ^A^,^'', can be calculated with a small parameter 77 as {4] 

Sh{u) = ^{h [p^, 4"^*] - h [po, Ko, 1^*0]) , (15) 
^AW(^) = ^ (a [Pv, - ^ [Po, Ko, K*]) , (16) 

^A(")(u;) = i (A [pt , 4-), 4+)*] - A [po, ^0, ^l]) ■ (17) 

Here, the FAM densities, prj and kI^\ are slightly changed from those at the ground state 
and calculated with the modified quasiparticle wave functions as 

p^ = {V* + r]UX){V + r]U*Yf 
4+) = {V* + riUX) {U + 7]V*YY (18) 
4"^ = iy* + r]UY*){U + 7]V*X*f. 



From these FAM formulae, we may calculate SH^u^ for a given vector {X,Y), without the 
explicit calculation of the comphcated residual interactions, v^^yyi. 

Now, we apply the same trick as we did in Sec. IIIBt to calculate the QRPA matrix 
elements. To obtain A^^^^v', we choose 

Xal3 = Sa^i'Spu', Yap = 0, (19) 

while B^y^^iyi is obtained by choosing 

Xal3 = 0. Yaji = Sa^i'SiSu'- (20) 

If we want {A± B), we may use the following 

Xa/3 = San'Spu'- Ya/3 = ±5a^_i' S j3u' , {A ± B) ^y^^i yi . (21) 

In this way, the QRPA matrix can be explicitly constructed by the FAM approach. 

IV. NUMERICAL RESULTS 

In this section, we show the numerical results of the QRPA based on the HFB ground 
state. The present FAM approach to the calculation of the QRPA matrix has been im- 
plemented in the computer code hfbrad for the spherical Hartree-Fock-Bogoliubov (HFB) 
calculation [l2]. The static HFB solution is obtained in the radial coordinate space dis- 
cretized with a mesh of 0.1 fm in the box of 20 fm. The maximum angular momenta for the 
neutrons and protons are set to be jmax = 21/2 and 15/2, respectively. We use the energy 
density functional of the SkM* parameter set and of the volume-type pairing 

V^pair(r - = Vo(l - P.)Kr- (22) 

with the strength Vq = —90 MeV fm'^. The quasiparticle energy cutoff is set at E^^ = 
60 MeV, unless otherwise specified. Then, the FAM construction of the QRPA matrix is 
performed. We include all the two-quasiparticle states in the quasiparticle space truncated 
by -Eqp. The number of two-quasiparticle states with the J'^ = 0+ is A'2qp = 1741 for 
= 60 MeV for ^^osn. 
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A. Monopole strength function 



We calculate the QRPA matrix using the FAM presented in Sec. Illli To obtain the QRPA 
matrix of A, we repeat the FAM calculation adopting vectors (X^^, Yq,^) = 0), with 

A^2qp different pairs of /iV. Then, to obtain the matrix B, we use another A^2qp kinds of 
vectors, (Xq,^,Yo;9) = {0,Safi'Si3u')- Once the A and B matrix is calculated, we resort to a 
routine ZGEEV in the Lapack libraries [13] to diagonalize the QRPA matrix in Eq. ([1]). The 
diagonalization of the QRPA matrix produces the normal-mode excitation energies Un and 
eigenvectors F*^")). This defines the n-th normal-mode creation operator, 

= E - • (23) 

fJ,>U 

Then, we calculate the transition matrix elements of a one-body operator F between the 
ground state |0) and the n-th excited state \n) = fin\0), as 

{n\F\0) = (HFB|[fi„,F]|HFB) = {X^;!.^* F'u + Yi^^* K") (24) 
and F^"^ are a^a"!"- and aa-parts of the one-body operator F defined by 

F = Y1 ^'^'^ = ^0 + 5^ + F;^X<^,) + J2 F>W (25) 

kl iJ,>u flu 

In case that the operator F is a hermitian operator conserving the particle number, they 
are given by 

Fll={U^FV*-V^F^U*)^^, (26) 

and = 

The result of the present approach to the isoscalar monopole strength function {F = 
'Ylik'^V) ^^°Sn is shown in Fig. [T] In this nucleus, the protons are in the normal phase, 
while the neutrons are in the paired superfluid phase with an average gap of 1.3 MeV. We also 
perform the iterative solution of the inear-response calculation (i-FAM) for a given energy 

following the procedure in Ref. |4[. We adopt an iterative algorithm of the generalized 

conjugate residual (GCR) method [ij]. The energy E is taken from to 45 MeV with 

^E = 0.1 MeV, to which we add an imaginary part 7/2, E — )■ E + i'~f/2. In this calculation, 

we adopt 7 = 1 MeV. To compare the results, we smear the transition amplitudes of Eq. ( 12^ 

in the m-FAM, with a Lorentzian function of 

m 7/2^ \{n\F\0)\' 

S{uj;F) = > — . (27) 

^ ^ TT ^ (u; - a;„)2 + 72/4 ^ ^ 
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FIG. 1: (Color online) Isoscalar monopole strength function for ^^^Sn calculated with m-FAM and 
i-FAM. The green dotted curve is obtained with the iterative method without explicit calculation 
of the QRPA matrix, while the red curve is the present approach to the QRPA matrix with a 
Lorentzian smearing of 7 = 1 MeV. 

It is clear from Fig. [T]that the two independent calculations give essentially the identical 
result. There is a small deviation in the peak strength near zero energy, which is due to the 
fact that, in the vicinity of the zero energy, the strength function obtained by the linear- 
response calculation with the complex energy E ^ E + i'j /2 is slightly different from the one 
using the smoothing with the Lorentzian form of Eq. ( 1271) . Actually, this peak is associated 
with the Nambu-Goldstone mode of the pairing rotation for neutrons. The energy of this 
spurious mode in ^^°Sn is calculated as cung = 662 keV. In principle, we should obtain 
i^NG = because the present calculation is fully self-consistent. However, the energy of the 
spurious mode is extremely sensitive to numerical errors in calculation of the QRPA matrix 
elements, which leads to the sizable energy shift. Note that other physical excitations are 
practically not affected by these errors. In fact, we have confirmed that this single state 
carries 99.9 % of the total strength associated with the neutron number operator, |(n|iV|0)p. 

In the i-FAM, we calculate the response function for the external field of the isoscalar 
monopole, however, we do not know the eigenenergies and eigenvectors of the QRPA normal 
modes. The present method provides us with this missing information. 
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B. Pair transfer modes in neutron-rich Pb isotopes 



A feature of the m-FAM different from the i-FAM is an exphcit calculation of the normal 
mode. This is useful for studies of low-lying excited states. In this respect, the m-FAM has 
an advantage over the i-FAM. In this section, we study the neutron-pair-transfer strength for 
low-lying excited 0^ states in Pb isotopes. The pair-transfer operators for the (L, S) = (0, 0) 
pair are defined as 

^ ^ /" dff{r)^i{fi)^^Ji{r\), (28) 



P = -— I rfr7(r)V>„(r t)^„(r ;), (29) 



'Atx 

for pair- addition and removal operators, respectively. Here, ip\{r,a) and '?/'„(r, a) indicate 
the creation and annihilation field operators, respectively, for neutrons with spin a at the 
position r. The radial form factor /(r) is chosen as unity in this paper. We calculate the 
transition probabilities from the ground state of the nucleus with the neutron number N to 
an J'^ = 0^ excited state \n) in the nucleus with ± 2, 

5(add;gs^n) = |(n|P^|0)p, (30) 
S(rem;gs^n) = |(n|P|0)|l (31) 

They are given by exactly the same expression as Eq. (12^ . For the pair addition mode, 
and P"^ are replaced by 

(^')?^ = ^ E - ui,ul) , (32) 

^ fc>0 

(PV^ = ^ E (^^/^^S- - ^^z^^^'^) ' (33) 

respectively. For the pair removal mode P, we have the same expressions by interchanging 
U ^ V with an opposite sign. 

We perform the calculation for even-even Pb isotopes from = 126 to A^ = 142. The 
HFB calculation with the SkM* parameters predicts the ground-state properties in Table [H 
The neutron pairing gap gradually increases as the neutron number. In Fig. |2l we show 
two-neutron-transfer strengths as a function of the (Q)RPA normal-mode excitation energy 
E. The ground state in ^°*Pb is in the normal phase with A„ = Ap = 0, due to its 
doubly closed-shell configuration. In this nucleus, the lowest mode around E = 2.6 MeV 
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TABLE I: Calculated chemical potentials A„ [ MeV ] and average pairing gaps A„ [ MeV ] for 
neutrons in Pb isotopes. The protons are in the normal phase, Ap = 0, for all isotopes. 



N 




An 


126 


-6.27 


0.00 


128 


-4.87 


0.77 


130 


-4.68 


1.03 


132 


-4.50 


1.19 


134 


-4.33 


1.31 


136 


-4.17 


1.42 


138 


-4.03 


1.53 


140 


-3.90 


1.62 


142 


-3.80 


1.70 



corresponds to the ground state in neighboring nuclei with ± 2. This corresponds to 
the pairing vibration which carries significant strengths both for two-neutron addition and 
removal modes; B{add]^^^Ph ^^^^ Pb) = 7.47 and 5(rem;208pb Pb) = 4.93. For the 
removal mode, the calculation suggests another 0"'' state with a sizable strength, located at 
about 1.4 MeV higher than the ground state in ^°^Pb. The excited 0"*" state at -Eex = 1.17 
MeV was observed by the {p, t) reaction whose cross section is about 10 % of that of the 



ground state [15 



lel. 

In the normal phase, the chemical potential A is not uniquely determined. Although the 
pair strengths are invariant, the energies of the normal modes {E^'^'^, E"^*^^) depend on the 
choice of the chemical potential in ^°^Pb. However, the sum of them, E^'^'^ + E^'^'^, does 
not depend on the chemical potential. In fact, the chemical potential for ^°^Pb in Table [T] 
is determined by the condition that the pair addition and removal modes have the same 
excitation energy. 

For superfiuid isotopes with the finite pairing gap (A„ 7^ 0), the energy in Fig. [2] can be 
regarded as 

E = Eex(iV ± 2) - E^,iN) t2X^ Eex(iV ± 2) - E^^iN ± 2). (34) 
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There are no significant pair-removal strength in excited O"*" states for ^^°~^^^Pb. The lowest 
state is calculated to be about 6 MeV higher than the ground state. In contrast, the pair- 
addition strength is present at low energy; E^,^ = 3.0 MeV for ^^°Pb and the energy is even 
lowered by increasing the neutron number, leading to the minimum value of E^^ = 2.3 MeV 
in ^^'^Pb. Further increasing the neutron number, this lowest 0+ transfer mode gradually 
changes from the pair-addition into the pair-removal character. In ^^^Pb, the lowest mode 
carries very little pair-addition strength. 



V. COMPARISON IN COMPUTATIONAL POINT OF VIEW 



A. Iterative FAM (i-FAM) 

The iterative algorithms adopted in the former i-FAM calculations 3|-l7| are slightly differ- 
ent from each other. Nevertheless, these calculations are composed of common ingredients: 

1. The (Q)RPA matrix is not explicitly constructed. 

2. The linear-response equation is solved for a fixed energy u = E + i'-f/2 with a given 
external field. 

Thus, the i-FAM solves a linear algebraic system of the form Mx = b. The algorithm jl^ 
involves matrix-by- vector multiplication, Mp for a given p, to produce a succession of vectors 
converging to the solution: 

(35) 

fn = b — MXn 0. 

The precision of the solution can be measured by e„ = and the iterative proce- 

dure stops when e„ < e. With this technique, we do not need to compute the values of the 
matrix elements of M , since it is enough to know the product of Mp for a given vector p. 
The GCR algorithm requires calculation of Mp twice at each iteration. From the view point 
of the memory resources, the method has a great advantage because the biggest arrays to 
be used in the code are of the order of 2A^2qp, instead of those for the matrix, 2A^2qp x 2A^2qp- 
For the present calculation for ^^"^Sn, we have employed the following parameters: the 
precision e = 10~^, the FAM parameter in Eq. ( |T8l) rj = 10~^, the energy range of < < 45 
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FIG. 2: Neutron pair transfer strengths to excited states in Pb isotopes. The energy is regarded as 
the approximate excitation energy in the nucleus with N ±2, except for ^°^Pb. The pair-addition 
strength is shown as the positive value, B(add;iVgs {N + 2)ex), while the removal strength is 
shown as the negative strength, — S(rem; A^gg {N — 2)ex)- These quantities are dimensionless, 
since the radial form factor is chosen as unity, /(r) = 1. Note that the following pair strengths 
are out of vertical range: 5(add) = 7.47 and B{vem) = 4.93 at E = 2.61 MeV in ^ospb and 
5(rem) = 3.74 at ^ = 5.95 MeV in ^lOpb. 13 



MeV discretized with a mesh of AE = 0.1 MeV with a smoothing parameter 7 = 1 MeV. 
The speed of the calculation depends on the number of iterations required to reach the 
convergence, which is particularly affected by the precision parameter. In the present case, 
the number of iterations ranges from about 30 at low energy, to about 160 where the strength 
has a peak (around 15-20 MeV). To start the iteration procedure, we need an initial vector 
xq. We used the solution obtained aX E — AE as the initial vector for the energy E. 



B. Matrix FAM (m-FAM) 

To construct the full QRPA matrix, we need to compute the right-hand side of Eq. ( IT3|) 
for 2A^2qp kinds of vectors {X,Y). In case that N2qp is relatively small, this is the most 
time-consuming part for the m-FAM. It also requires the memory capacity to store the 
matrix of order of 2A^2qp x 2A^2qp- However, there is a practical advantage in the m-FAM; 
the calculation of the residual fields, SH'^'^, is easier in the m-FAM than the i-FAM. For a 
given vector (X, F), we calculate the densities, {pj^^ni^^), in Eq. flTSl) . In the i-FAM, since 
the vector (X, Y) is updated every iteration, we construct by the matrix operation as 

(Pnh = E + ^ E (^^^- + ^ E ^7^- j ' (36) 

and we have similar expressions for k'^^^ Here, we need to sum over the three quasiparticle 
indices (a,/3,7). In the case of m-FAM, we can omit these summations, because of the 
simple form of {X,Y) such as Eqs. f|T9l) and fl20l) . For instance, the calculation of A^^^fj_'^i 
can be carried out by adopting Eq. (fTOl) . Then, the densities are simply given by 

(4^^)^, = /^o + r?f/vf^.v, (37) 

Thus, the calculation of these densities is faster in the m-FAM than in the i-FAM. In the 
present numerical calculation, we have found that this reduces the computation time for 
building the QRPA matrix, becoming 1/5 of the original time with Eq. fl5Bl) . 

After the calculation of the QRPA matrix, we diagonalize the matrix using the Lapack 
libraries. The computational task required for this diagonalization scales as (2A^2qp)'^, while 
that for the calculation of the QRPA matrix elements scales as as (2iV2qp)^. Thus, increasing 
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TABLE II: Computational time of the iterative FAM (i-FAM) for the monopole strength function 
and the matrix FAM (m-FAM) for all the J'" = 0+ eigenstates with different cutoff energies (-Eqp), 
relative to that of i-FAM with E^^ = 60 MeV. The size of the QRPA matrix is determined by the 
number of 0"*" two-quasiparticle states (A''2qp) in the model space, which is varied by changing E^. 





2 X iVsqp 


i-FAM 


m-FAM 


60 MeV 


3482 


1 


0.16 


80 MeV 


4656 


1.43 


0.38 


100 MeV 


5842 


1.93 


0.60 


120 MeV 


7156 


2.64 


1.26 


140 MeV 


8336 


3.27 


1.77 


160 MeV 


9528 


4.08 


2.56 



the quasiparticle model space, the diagonahzation will eventually become a hot spot in the 
computation. 

C. Comparison of computational time for ^^'^Sn 

In Table [Tll we show the relative CPU time of the calculation of monopole strength in i- 
FAM and the QRPA calculation with the m-FAM, increasing the number of two-quasiparticle 
states (increasing -Eqp). In the present calculation with the quasiparticle-energy cutoff of 
-Eqp = 60 MeV, the m-FAM is six times faster than the i-FAM. The computation time for 
i-FAM shows a weak scaling with respect to the size of the model space, close to the linear 
dependence on A^2qp- On the other hand, the computation time for m-FAM indicates a 
scaling between A^lqp ^2qp- Therefore, for treating the larger model space, the i-FAM 
has a computational advantage over the m-FAM, with respect to this scaling property and 
the memory requirement. 
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D. Parallelization 



In contrast to the present FAM calculations for spherical nuclei, the FAM calculation for 
deformed systems may require significant computational resources with modern massively 
parallel supercomputers. In the case of the i-FAM, the calculation of the strength functions 
at a given energy is independent from the other energies. This leads to an obvious kind 
of parallelization with the number of processors equal to that of energy points. Typically, 
the number of energy points is order of 100 at most. Since the iteration process cannot be 
parallerized, the parallel computation with processors beyond this number is not trivial. 

In contrast, the large-scale parallel computation can be easily done in m-FAM. The m- 
FAM calculates the matrices, A^y^^iyi and B^y^^iyi. For a given pair of yuV, the matrix 
elements of A^^y^i are produced from the vector of Eq. f|T9|) . and those of B^yy^i are from 
that of Eq. ( 120|) . Since these procedures are completely independent, it is easy to utilize 
the processors whose number is as large as 2 x A^2qp- For deformed systems, the number of 
two-quasiparticle states could be order of 10^, even assuming the axial symmetry [ill, 
Thus, the massive parallelization of this magnitude can be achieved in the m-FAM. 



VI. CONCLUSIONS 



We have presented a feasible method to construct the matrix of the RPA and the QRPA, 
based on the idea of the finite amplitude method (FAM). Since all the residual interactions 
are numerically estimated by the FAM, it does not require complicated programming. The 
method can be easily implemented with existing HF and HFB codes, to turn them into the 
RPA and QRPA codes. We have used the hfbrad code to test the present method. An 
advantage of the present method, called matrix FAM (m-FAM), over the i-FAM is that the 
m-FAM provides the normal-mode vectors explicitly. 

In the computational aspect for a large-scale problem (increasing A'"2qp), a disadvantage 
of the m-FAM over the i-FAM is that the computational task scales as (A^2qp)^ for the 
calculation of the matrix elements and as (A^2qp)^ for the diagonalization. In addition, it 
requires a large memory capacity to store the QRPA matrix of 2N2qp x 2iV2qp. On the 
other hand, there are some advantages as well. For problems of small dimensions, the 
computational time is shorter than the former iterative FAM (i-FAM). For problems of 
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large dimensions, the m-FAM may resort to a massively parallel computer, because the 
calculation of the matrix elements can be easily parallelized up to the number of processors 
equal to twice of the number of two-quasiparticle states. Thus, the m-FAM could be a strong 
candidate of the new QRPA code for deformed superfiuid nuclei, for the use of the massively 
parallel computers. 

In many respects, the i-FAM and the m-FAM have complimentary character to each other. 
In addition, there are intriguing developments in the iterative approaches to the lo w-ly ing 



normal-mode solutions: The conjugate gradient method was adopted in Refs. [18|, ll9| for 
the RPA solutions in the 3D real-space representation. Very recently, the iterative Arnoldi 
method, developed in Ref. j2o|, has been applied to the calculation of the lowest solutions 



of the QRPA in spherical nuclei 



21| . These novel technologies and the present FAM may 



provide a powerful tool for developments of the efficient and feasible (Q)RPA codes for 
deformed systems. 
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